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For plane-wave and many-spiral states of the experi- 
mentally based Luo-Rudy 1 model of heart tissue in large 
(8 cm square) domains, we show that an explicit space-time- 
adaptive time-integration algorithm can achieve an order of 
magnitude reduction in computational effort and memory — 
but without a reduction in accuracy — when compared to an 
algorithm using a uniform space-time mesh at the finest res- 
olution. Our results indicate that such an explicit algo- 
rithm can be extended straightforwardly to simulate quanti- 
tatively large-scale three-dimensional electrical dynamics over 
the whole human heart. 
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Understanding the dynamics of excitable media such 
as heart tissue is a problem of substantial interest to 
physicists, physiologists, biomedical engineers, and doc- 
tors. For reasons not yet understood experimentally, 
the healthy time-periodic spatially-coherent beating of 
a human heart will sometimes change to a nonperiodic 
spatially-incoherent fibrillating state in which the heart 
cannot pump blood effectively (leading to death if suit- 
able treatment is not administered quickly). It would be 
valuable to understand how the onset of arrhythmias that 
lead to fibrillation depends on details such as the heart's 
size W, geometry, electrical state, anisotropic fiber struc- 
ture , and inhomogeneities. A deeper understanding of 
the heart's dynamics may also make possible the inven- 
tion of protocols by which electrical feedback could be 
used to prevent fibrillation ||^. 

Because of many experimental difficulties in studying 
the three-dimensional dynamics of a heart simula- 
tions of cardiac tissue (and more generally of excitable 
media) play an extremely important role in identifying 
and testing specific mechanisms of arrhythmia. However, 
quantitatively accurate simulations of an entire three- 
dimensional human heart are not yet feasible. The essen- 
tial difficulty is that human heart muscle is a strongly ex- 
citable medium whose electrical dynamics involve rapidly 
varying, highly localized fronts (see Figs. |l| and I). The 
width of such a front is about 0.05 cm and a simula- 
tion that approximates well the dynamics of such a front 
requires a spatial resolution at least 5 times smaller, 
Aa; « 0.01 cm. The muscle in an adult human heart has a 
volume of about 250 cm^ and so a uniform spatial resolu- 
tion of 0.01 cm would require a computational grid with 
w 3 X 10^ nodes. Depending on the assumed material 



properties of the heart and on the quantities of interest to 
analyze, up to 50 floating point numbers might be associ- 
ated with each node, requiring the storage and processing 
of about 10^" numbers per time step. The fastest time 
scale in heart dynamics is associated with the rapid de- 
polarization of the cell membrane, about 0.1 msec in du- 
ration, and a reasonable resolution of this depolarization 
requires a time step about a fifth of this. At « 0.02 msec. 
Since arrhythmias such as fibrillation may require several 
seconds to become established, the 10^*^ numbers associ- 
ated with the spatial mesh would have to be evolved over 
about 10^ time steps. Such a uniform mesh calculation 
currently exceeds existing computational resources and 
has not yet been carried out. 

A clue about how to improve heart simulations comes 
from experiments and simulations which suggest 
that the electrical membrane potential V{t, x) in the fib- 
rillating state consists of many spirals (for approximately 
two-dimensional tissue such as the atrium, see Fig. ||) or 
of many scroll waves (for thicker cardiac tissue such as the 
left ventricle |^ ) . A striking feature of these spatiotem- 
poral disordered states is that the dynamics is sparse: at 
any given time, only a small volume fraction of the ex- 
citable medium is occupied by the fronts, and away from 
the fronts the dynamics is slowly varying in space and 
time. It may then be the case that the computational 
effort and storage can be greatly reduced, from being 
proportional to the volume of the excitable medium (the 
case for a spatially uniform mesh) to being proportional 
to the arclength (in 2d) or surface area (in 3d) of the 
fronts. 

In this Letter, we show for representative solutions of 
the quantitatively accurate Luo-Rudy 1 (LRl) membrane 
model of cardiac tissue [^) that an explicit space-time 
adaptive-mesh-refinement algorithm (AMRA) [Q can in- 
deed take advantage of the sparse excitable dynamics to 
reduce by an order of magnitude the computational ef- 
fort and memory needed to simulate arrhythmias in large 
domains. Further, we show that there is no significant 
reduction in accuracy when using an AMRA compared 
to an algorithm that uses a spatially uniform mesh at 
the finest resolution of the AMRA. Since the AMRA is 
explicit in time and has a fairly simple data structure 
consisting of nested patches of uniform Cartesian meshes, 
the AMRA can be parallelized straightforwardly |^ , lead- 
ing to a further reduction in computational effort by the 
number of processors. The AMRA is also general and 
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does not require the use of reduced models [^[Q], which 
increase efficiency but sacrifice experimental accuracy by 
using fewer variables and perhaps explicitly eliminating 
rapid variables. The results presented below suggest that 
a quantitatively accurate AMRA simulation of fibrilla- 
tion in an entire human left ventricle for several seconds 
with an effective 0.01 cm resolution should already be 
practical with existing computers. This is highly encour- 
aging since further improvements to such algorithms are 
possible as discussed below. 

In the following, we discuss some details of the AMRA 
and then its accuracy and efficiency for simulations of the 
LRl model in large one- and two-dimensional domains. 
Our particular algorithm was a straightforward modifica- 
tion of an AMRA that has been used by other researchers 
to integrate hyperbolic sets of partial differential equa- 
tions such as the Euler equations of fluid dynamics . 
Since key mathematical and algorithmic details are avail- 
able elsewhere Q , only some essential ingredients and our 
modifications of them are briefly described here; a more 
detailed discussion will be given elsewhere . 

The AMRA approximates a given continuous field such 
as the cardiac membrane potential V{t,x.) on a set of 
nested locally-uniform patches of d-dimensional Carte- 
sian meshes in a d-dimensional Cartesian box On 
each patch, spatial derivatives in the dynamical equations 
are approximated by second-order-accurate finite differ- 



ences and an explicit method (we use forward-Euler |10|) 
is used to advance in time. The power of the algorithm 
arises from its ability to automatically and efficiently re- 
fine or coarsen the representations of fields by varying 
the number of grid points locally to achieve a specified 
truncation error. A further reduction in computational 
effort is achieved by allowing the the time step to change 
locally with the spatial mesh 0. In related prior work, 
Quan et. al. ||l^ have studied cardiac models using spa- 
tially adaptive time steps but with a uniform spatial mesh 
and alternation of implicit and explicit time steps, while 
Moore ||l^ has studied reaction-diffusion equations us- 
ing a spatially-adaptive fully-implicit method but with a 
spatially-uniform adaptive time step. To our knowledge, 
ours is the first study of an algorithm for excitable me- 
dia for which the spatial and temporal resolutions change 
locally. 

An important subtlety is that our AMRA was designed 
for hyperbolic equations but is here applied to an ex- 
citable medium which is described by parabolic equations. 
For explicit time integrations of hyperbolic equations, the 
Courant-Friedrichs-Lewy (CFL) condition for the onset 
of numerical instability bounds the largest possible 
local time step At by the first power of the local spa- 
tial resolution Ax. For parabolic equations, the stability 
condition for an explicit algorithm bounds the time step 
by Ax^ for Aa; sufficiently small. In the LRl model, 
this dependence is evident only for spatial resolutions an 
order of magnitude finer {Ax < .0025 cm) than those 



used in our calculations. For resolutions in our range 
of interest, the fast reaction kinetics, not the diffusion 
operator, sets the stability limit on the time step 
A standard way to avoid the stability restriction on At 
is to use a semi- or fully-implicit time-integration algo- 
rithm [pl prlJl^ . We have estimated that by using an ex- 
plicit integration scheme, our time steps on the finest 
meshes are about an order of magnitude smaller than 
those needed to achieve a 10% relative error in the speed 
of the front (AMRA uses 0.003 ms as opposed to the 
value 0.04 ms for the semi-implicit case) Q. However, 
one cannot conclude that a semi-implicit algorithm is 
automatically better than our explicit one since, for a 
fixed spatial resolution, the larger time step allowed by 
a semi-implicit method may give less accuracy during 
the upstroke |13| and require more computation (some 
of these issues will be discussed quantitatively elsewhere 
for the Id case [^). Since the spatiotemporal dynam- 
ics of even the most detailed cardiac membrane models 
are not yet understood and the relation between speci- 
fied local truncation error and correct dynamics is also 
not understood, the present calculations should be con- 
sidered as an early but significant step in finding a good 
balance between efficiency and accuracy for simulating 
arrhythmias in large domains and over long times. 

Our results for the AMRA were obtained for the quan- 
titatively accurate LRl model |^, which in 2d can be 
written in the form: 

CmdtV{t, x,y) = ^ {gxdlV + gydyV) - /ion(m) - Istim{t, x, y), 



where F(i,x) is the membrane potential at time t and 
at position x = {x,y), Cm is the membrane capacitance 
per unit area, /3 is a surface-to-volume ratio of a heart 
cell, gx and gy are membrane conductivities (generally 
not equal since the heart is anisotropic), /ion is the total 
ionic current flowing across the membrane, and /stim is a 
specified current injected to initiate a propagating wave. 
(For all calculations reported below, the boundary con- 
dition (n ■ y)V = was used, where fi is the unit vector 
normal to a given boundary point.) The seven voltage- 
sensitive membrane variables TOi(t,x) for the LRl model 
determine the fiow of various ions across the membrane 
and satisfy ordinary differential equations, which are also 
integrated by a forward-Eulcr method. The same mem- 
brane parameter values as those of Rcf. |^ were used 
except for the calcium conductivity gca. in the /jon term, 
whose value was changed from 0.09 to 0.045 (in units of 
mf2~^ • cm~^). The medium was isotropic with g^ and 
gy set to 1 kn~-^ ■ cm~^ and f3 set to 3000 cm~^. These 
values shortened the action potential duration and led 
to dynamical states with many spirals, providing a more 
challenging test of the AMRA. 

In addition to the physical parameters in Eq. (pi), many 
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numerical and algorithmic parameters need to be speci- 
fied 0;^. Several of the more important choices are an 
initial resolution for a uniform coarse mesh covering the 
domain (we used Aa; = 0.05 cm), the temporal resolution 
for the coarse mesh (we used At = 0.012ms), the max- 
imum number of grid levels allowed for refinement (we 
used the value 3), the factor by which the spatial mesh 
is refined locally (we chose the factor 2), the error tol- 
erance used in the Richardson extrapolation estimate of 
the local truncation error (we chose e = 2 x 10~^); and 
the number of time steps to elapse before estimating a 
local error and regridding (we chose 2). 

As a first demonstration of the effectiveness of the 
AMRA, Fig. summarizes a 3-level calculation of the 
LRl model in a Id domain of length L = 9 cm. The sys- 
tem was stimulated at t = with a 0.2 cm square pulse 
along the left edge of the domain, which evolved into a 
front propagating to the right (the spatial profile is in- 
dependent of the initial condition and of the system size 
for L > 9 cm) . One can see from the spatial profile in 
Fig. ^ at time t = 240 ms how narrow is the front (re- 
gion of depolarization) compared to the profile's extent 
and this specifically is what makes numerical simulation 
of highly excitable media so difficult. In the vicinity of 
the front, Fig. shows the grid structure which was au- 
tomatically calculated by the ARMA; the colors black, 
green, and red indicate the coarse, fine, and finest mesh 
regions respectively. Taking into account the reduction 
of spatial mesh points and the asynchronous updating 
of grid points using spatially varying time steps , the 
AMRA overall used a factor of 3.6 fewer grid points and 
did less computational work by a factor of 9 for the LRl 
model than a constant-time-step uniform-spatial-mesh 
forward-Euler code using the finest space-time resolu- 
tions of the AMRA. The spatial adaptivity of the time 
step accounts for a factor of 2 in this factor of 10 and 
so is an important part of the algorithm. The temporal 
profiles at a fixed point in space, the front speeds, and 
the times between peak and recovery at a fixed point in 
space (action potential duration) for the AMRA and for a 
high- resolution uniform- mesh code (discussed in Ref . j|] ) 
agree within 0.1% relative errors except at the peaks of 
the temporal profiles, where the relative error is about 
4%. We conclude that there is no significant loss of ac- 
curacy when using the more efficient AMRA. 

Fig. H shows how the AMRA performs for the 
LRl model in a large square domain of size L = 8 cm, 
using the same parameter values as the Id case, for which 
spirals are unstable and break up into other spirals. This 
complex many-spiral dynamical state is a much stronger 
test of the efficiency and utility of an AMRA than Fig. |l| 
since the geometry of the fronts fluctuates strongly in 
time. A multi-spiral state was initiated by a standard 
S1-S2 stimulation protocol |5j in which a right-going pla- 
nar pulse is created by stimulating the left edge of the do- 
main (the SI stimulus), and the lower left quadrant of the 



domain is excited (the S2 stimulus) 334 ms later, when 
the left half of the domain has returned to rest but the 
right half is still repolarizing. A comparison of the field V 
with the instantaneous grid structure approximating V 
is given in Fig. |^ 1346 ms after S2 and demonstrates how 
the AMRA is able to increase automatically the space- 
time resolution only in the vicinity of the fronts, greatly 
decreasing the overall computational effort since, at any 
given time, the sharp fronts indeed occupy only a small 
fraction of the domain. The total number of mesh points 
used by the AMR varies substantially with time, from 
3 X 10* to 7 X 10* mesh points with an average of 5 x 10*. 
A comparison of these results with those required by a 
uniform-spatial-mesh constant-time-step code using the 
finest AMRA resolution [H shows that the AMRA uses 
about 8 times fewer mesh points, requires less integration 
work by a factor of 12, and achieves a speedup of about 
a factor of 11 

The above results can be used to estimate the com- 
puter time needed by the ARMA to integrate for one 
second the LRl model for a 3d section of left ventricular 
wall of dimensions 8 cm x 8 cm x 1 cm, with an effective 
fine uniform mesh resolution of Ax = 0.0125 cm in space 
and At = 0.003 msec in time. On a Pentium III 500 MHz 
computer, we found that a 3-level 2d AMRA calculation 
at this resolution takes about 3 days. The time for the 
3d calculation then can be estimated by assuming that 
each of the spirals in Fig. ^ becomes a continuous stack 
of spirals (a scroll wave) , with the stack transverse to the 
square sides of the domain , and correspondingly that 
the mesh refinements extend uniformly from the 2d case 
through the transverse direction. A 3d AMRA calcula- 
tion should then take roughly 15 days, which is a factor 
of 17 speedup over the 9 months required to complete 
a similar calculation using a uniform space-time mesh 
with the above resolution. Without substantial change 
to the AMRA, an additional speedup of at least 10 can 
be gained by using a distributed parallel computer with 
100 Pentium III processors, and another speedup of 5 
by using table-lookups to avoid the many exponentia- 
tions associated with the integration of the membrane 
variables mi[t). These further gains would reduce the 
total simulation time for one second of the LRl model in 
this 3d domain to 7 hours or less. (With a substantial 
modification to make the AMRA semi-implicit, another 
reduction by a factor of 2-3 might be possible.) Simula- 
tion of an entire heart (a factor of 4 greater in volume) 
for one second with a LR 1 model should then be pos- 
sible on the time scale of one day, which is acceptably 
fast for exploring many interesting questions about the 
dependence of arrhythmias on parameters. 

In summary, we have shown that an explicit space- 
time adaptive algorithm |0] using one of the simplest pos- 
sible data structures (a hierarchy of Cartesian meshes) 
can already attain an order of magnitude reduction in 
computational effort and memory when applied to the 
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experimentally based LRl cardiac membrane model 
and that this reduction is achieved without incurring a 
corresponding reduction in accuracy when compared to 
an explicit code using a uniform space-time mesh. Im- 
portant next steps include determining whether the algo- 
rithm can be improved by using implicit time integration, 
generalizing the method to curved boundaries, and mak- 
ing specific applications to the initiation and control of 
human arrhythmias. 

We thank M. Berger, Z. Qu, and A. Garfinkel for useful 
discussions and especially M. Berger for making available 
to us one of her AMRA codes. This work was supported 
by a NSF Graduate Research Fellowship, by NSF grant 
DMS-9722814, and by NIH grant R29-HL-57478. 
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FIG. 1. (a) Spatial profile V{t,x) at time t = 240 ms 
for a Id front propagating to the right in a domain of 
length L — 9 cm, as calculated by a 3-level adaptive mesh 
refinement algorithm (AMRA) for the Luo-Rudy 1 (LRl) 
cardiac model [^. The three regions of coarse, fine, and 
finest mesh resolution (from Ax — 0.05 cm. At = 0.012 ms 
to Aa; = 0.0125 cm, At = 0.003 ms) are indicated by the 
black, green, and red portions of the curve, (b) Blowup of 
the small interval indicated near x — 8.4cm in (a), showing 
how the 3-level mesh structure (vertical lines) has automati- 
cally resolved the sharp front. 



FIG. 2. (a): Three-level AMRA calculation of the 2d 
LRl model at time t = 1346 ms after stimulus S2, in a 
square domain of length L = 8 cm. Field value ranges 
for V{t,x,y) are color coded with blue for V > — 5mV, red 
for —5<V< — 65mV, and yellow for V < — 65mV. Param- 
eter values are the same as in Fig. |^. (b): The hierarchical 
Cartesian meshes of the AMR algorithm corresponding to the 
snapshot of V in (a). The yellow and green regions corre- 
spond to the fine (level 2) and finest (level 3) grids and track 
closely the fronts. 
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